A quantum Monte-Carlo method for fermions, free of discretization errors. 
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In this work we present a novel quantum Monte-Carlo method for fermions, based on an exact 
' decomposition of the Boltzmann operator exp(~f3H). It can be seen as a synthesis of several 

related methods. It has the advantage that it is free of discretization errors, and applicable to 
general interactions, both for ground-state and finite-temperature calculations. The decomposition 
is based on low-rank matrices, which allows faster calculations. As an illustration, the method is 
applied to an analytically solvable model (pairing in a degenerate shell) and to the Hubbard model. 
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Quantum Monte-Carlo methods offer an interesting way to obtain numerical results for large quantum systems 
|r-||. A number of Monte-Carlo methods, that go by names as auxiliary- field shell- model grand-canonical 
|3| or projection quantum Monte-Carlo |^ are based on the decomposition of the Boltzmann operator exp{—[3H) 
as a sum or integral over exponentials of one-body operators. The latter are easy to handle numerically. Simple 
\ algebraic expressions exist to calculate their grand-canonical or canonical trace |5|Q| or their overlap between 
Slater determinants |^,^. The sum or integral is then evaluated using Monte-Carlo techniques, most often Markov- 
chain Monte-Carlo techniques such as the Metropolis algorithm [||. The basic ingredients of such a decomposition 
are the Suzuki- Trotter decomposition pO| ] to separate non-commuting parts of the Hamiltonian and the Hubbard- 
' Stratonovich transform ||ll|,|l^ to linearize the two-body part of the Hamiltonian. Both ingredients lead to systematic 
(— I discretization errors in the calculations. Furthermore, for general two-body interactions these methods require many 
Q manipulations with dense matrices and hence a lot of CPU time. 

^ O ^ A finite-temperature method that is free of discretization errors was presented by Cerf jlSj. It is based on the 

expression 
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^Nj However, this method is limited to very specific types of interactions. It has only been applied to a model with 
pairing in a degenerate shell. Cerf 's method is not based on a decomposition into a sum of exponentials of one-body 
' operators. Hence the simple algebraic expressions for the (grand-)canonical trace and for the overlap between Slater 
determinants cannot be applied here. 

In this letter we present a new decomposition of the Boltzmann operator that is free of discretization errors, just 
like Cerf 's method, but that is applicable to general interactions and moreover results in a (partially continuous) sum 
over exponentials of one-body operators, just like auxiliary-field quantum Monte-Carlo methods. It can be seen as a 
^ , synthesis of the two. Furthermore this method is based on low-rank matrices, which results in a considerable speed-up 
'"^ in comparison with standard auxiliary-field quantum Monte-Carlo methods. The decomposition is an exact one, so 
the only error in the calculations is the statistical error originating from the Monte-Carlo sampling of the terms in 
O the decomposition (apart from the round-off error due to the limited machine precision). Because the method is 
^ ' applicable to general fermionic two-body interactions, it can be of use for nuclear, atomic as well as condensed-matter 
^ physics problems. However, it is not free of sign problems at low temperatures. But just like the standard Quantum 
Monte-Carlo methods, there is no sign problem for the attractive Hubbard model, the repulsive Hubbard model at 
rS , half filling |^ and the mean-field plus pairing model for even-even atomic nuclei jlBl-p^ . 

' The trick to arrive at the decomposition of the Boltzmann operator, is to replace —jSV in expression |l] with fi — fiV , 
where /.t is an arbitrary, real positive parameter. Adding to the exponent has no influence on the properties calculated 
with the Boltzmann operator, and it can simply be corrected by a factor . Instead of expression |] we now obtain: 
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This last expression is reminiscent of the expression for the partition function in the interaction representation derived 
in [0 , with the difference that here the factors have the form 1 — j^V{tp) instead of —V{t) (where V{t) is the two-body 
Hamiltonian in the interaction representation). It is this extra constant 1 that allows to make a decomposition into 
a sum of exponentials of one-body operators. To achieve this we start by constructing a decomposition for 1 — ^V. 
Hereby we build on expressions derived in reference | |l7|] . 
For the pairing Hamiltonian one has 



V 



-G 



(3) 



k,k'>0 



where the operator creates a particle in the corresponding single-particle state and with k the time-reversed state 
of the state k. The notation k,k' > denotes that the summation for k and k' should run over states with angular 
momentum projection > only. Using lemma 1 from reference we obtain 
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with 7 = — and Q half the number of single-particle states. Ak is the row matrix with a one on the entry 
corresponding to the state fc, and zero's anywhere else. As defined in reference p^ , for a square matrix Q, the operator 
0{Q) transforms a Slater determinant '^m, represented by the matrix M, into the Slater '^m', with AI' = Q M. If Q 
is non-singular, then 0{Q) is the exponential of a one-body operator. Note that this decomposition has a symmetry 
between the states fc > and their time-reversed states. This symmetry prevents sign problems for even particle 
numbers. 

For the repulsive Hubbard model one can take 
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Then 
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provided that cosh(7) = 1 - 
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with Ns the number of lattice sites and h^i 
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Lemma 1 from reference ll'i 



allows one to construct analogous exact decompositions of 1 — j^V ^ based on matrices of low rank, for any fermionic 
interaction. 

The Monte-Carlo algorithm has to sample over all values of m between and infinity, over all possible sets < 

h 



ii < • • • < im < 1 and at each interval over all the terms in the decomposition oi 1 — j^V . To perform this sampling 
numerically, a large number of intervals is taken. Let be the number of intervals. To each interval i we assign a 
fraction t; of the inverse temperature, such that X^il^i ti — I, and an index li to indicate that part of the decomposition 
of 1 — jjV that is inserted at that interval. If no part is inserted, then /j = 0. In total there are m out of intervals 
with li 7^ 0. Thus a term of order m in the decomposition ^ is represented by a configuration with m intervals for 
which li ^ 0. The sum of the coefficients between the j*'^ and the (j -|- 1)*'' insertion of 1 — -V , has to be equal to 



ij+i — tj. This scheme is visualized in figure |l|. 
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FIG. 1. A schematic way to represent the terms in the decomposition ^ The array Ti represents the inverse-temperature 



intervals, the array li represents the parts oi 1 — j^V that are inserted. 
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This representation is not unique. To obtain every combination with the right weight, we have to take into account 
an extra weight factor {N^ — m)\/Nx\ f^or a configuration of order m. The operator corresponding to a particular 
configuration can then be calculated as 



i=l 
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with Qj the j*'* part in the decomposition of 1 — ^V^, and Qo = 1- From a computational point of view it is 
advantageous to work in the interaction representation, i.e. 

Ul,r={^QlAU)^e-P^\ (8) 

with ti = ^''^d Qj{t) the j"* part in the decomposition of 1 — ^V[tl3). This decomposition is obtained by 

multiplying the row matrices used in the decomposition of 1 — —V, with the matrix e^*'^^°, where Hq is the matrix 

representation of the one-body operator Hq in the single-particle space. The operator Ui^r is a product of exponentials 
of one-body operators and thus the exponential of a one-body operator itself. Therefore, one can easily calculate its 
grand-canonical H or canonical trace [Q, or apply it to a Slater determinant. In this way, exact variants of the 
grand-canonical shell model § and projector quantum Monte-Carlo method [HJ] are obtained. When applying 
the grand-canonical or projector variant, fast updating techniques analogous to the ones presented by White et al. 
^ can be applied. The rank-two structure in the decomposition of the factors 1 — j^V{tf3), allows to make a quick 
update requiring only 8N^ flops, even for general interactions. These updates have to be performed only when li ^ 0, 
On average, this amounts to {m)MC times to update the whole configuration, where {m)MC is value of the order m 
in the decomposition ||, averaged over all Monte-Carlo samples. For the canonical algorithm such a fast updating 
technique is not possible. There the performance can be drastically improved using guided sampling p8[ |. Instead 

of the canonical (iV-particle) weight wi^r = tJ-"^{Nx — m)!Tr7v (^^,7-) /^x^-j one uses a local approximation lii^r that 
allows fast updates. After a number of steps these updates are then accepted or rejected collectively, according to 
the ratio q = ^p-^^jjj^- Using a generalized Metropolis algorithm jl^, that includes the factor /i"' '•^^™''' in the 

proposition probability, one can set up a very efficient Markov chain, with acceptance rates close to unity and with 
autocorrelation lengths of a few sweeps. 

Because the updating procedure is the most time-consuming step of the algorithm, the required CPU-time will 
be proportional to {m)MC- Therefore it is important to have a good estimate of this quantity in advance. If the 
grand-canonical or canonical algorithm is used, and if the weight wj^t is positive for all configurations, then one can 
show that 

{m)MC=^^-p{V)p■ (9) 

This shows that the CPU-time is proportional to the parameter fi, but also to the thermal expectation value {V)p of 
the residual interaction. Though fi is an arbitrary parameter, we have experienced that a good balance between low 
CPU-cost per sweep and fast mixing of the Markov chain is obtained by taking /i ~ (3\{V)i3\. Note that one has to 
take such that it is always larger than the largest value of m encounterd during the Monte-Carlo sampling. Apart 
from that, the method and the CPU-time it requires is independent of N^- Expression || can also be used to obtain 
a value for {V)i3 from the Monte-Carlo calculation. If wj,t can become negative, then one has to take into account 
the sign of w/,r in expression ^. To calculate expectation values for other observables, one has to use the techniques 
developed for auxiliary-field quantum Monte-Carlo methods, as described in references p|,|8|,p^ . 

To test our quantum Monte-Carlo method, we have applied it to a model with pairing in a degenerate shell. This 
many-body problem can be solved analytically using the seniority scheme pO| . We considered a system with 10 
particles in a degenerate shell {Hq = 0) of 20 single-particle states. We took V as in expression with G = IMeV. 
Figure || shows the results for the energy and the specific heat of the model, in the canonical ensemble. They agree 
perfectly with the analytical results. These observables were evaluated using equations 11 and 14 from reference [ p5[ . 

To demonstrate the versatility of the method, we have also applied it to the one-band repulsive Hubbard model. 
We have calculated the energy, specific heat and Coulomb energy in the canonical ensemble with 28 spin-up and 28 
spin-down particles, for a system with nearest- neighbour hopping only (strength t = 1), on an 8-by-8 lattice with 
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periodic boundary conditions and an interaction strength U = 2\t\. The results are shown in figure |3[ The average 
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sign of the weight wj^r for this system is shown as a function of the inverse temperature f3 in figure ^. All calculations 
were done on a Pentium Pro 200 MHz PC and a Digital Alphastation 255/300MHz workstation. 

Just like other quantum Monte-Carlo methods for fermions, our method is not free of sign problems at low tem- 
peratures. However, for a number of systems, a symmetry exists that guarantees good sign properties. This is among 
others the case for the attractive Hubbard model, the repulsive Hubbard model at half filling and the mean-field plus 
pairing model for even-even atomic nuclei. For other systems, one can expect the sign properties to be different from 
standard quantum Monte-Carlo methods, because a different decomposition for the Boltzmann operator is used. This 
might or might not be an improvement, depending on the specific situation, but has to be studied more deeply. 

The method can be improved by performing a (finite temperature) Hartree-Fock calculation for the Hamiltonian 
H in advance. This results in a one-body part Hq and a residual interaction V such that |(V)/3| is reduced. From 
expression ^ it is clear that this reduces {m)MC and hence reduces the CPU-time. Furthermore, one can expect that 
this will improve the average sign. 

In conclusion, we have developed a quantum Monte-Carlo method for fermions based on an exact decomposition 
of the Boltzmann operator into a continuous sum over exponentials of one-body operators. Because of the low-rank 
matrices used in the decomposition, fast matrix multiplications and efficient updating procedures can be applied. The 
method can be applied to systems with general two-body interactions, in the grand-canonical or canonical ensemble, 
or with ground-state projection. It allows to calculate thermal and ground-state expectation values for any observable. 
It is an exact method, apart from statistical errors. 

The authors are grateful to the F.W.O. (Fund for Scientific Research) - Flanders and to the Research Board of the 
University of Gent for financial support. 
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(a) 



FIG. 2. Energy E and specific heat C as a function of tem- 
perature T, for a model witli pairing in a degenerate sliell, 
as described in tlie text. Tfie curves correspond to tlie ana- 
lytical results, the errorbars indicate 2 x 2a intervals for the 
Monte-Carlo results. 
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FIG. 3. Total energy E and Coulomb energy 

Ec ~ {V)i3 + UN/2, (a), and specific heat C, (b), as a function 
of inverse temperature f3, for the Hubbard model described in 
the text, with 28 spin- up and 28 spin-down particles on an 
8-by-8 lattice, U — 2\t\. Units were chosen such that t — 1. 
(2(7 errorbars) 
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FIG. 4. Average sign {s}mc as a function of inverse tem- 
perature /3, for the same model as figure ^. {2a errorbars) 
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